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Abstract 

A new technique for calculating the time-evolution, correlations and 
steady state spectra for nonlinear stochastic differential equations is pre- 
sented. To illustrate the method, we consider examples involving cubic 
nonlinearities in an N-dimensional phase-space. These serve as a useful 
paradigm for describing critical point phase transitions in numerous equi- 
librium and non-equilibrium systems. The technique presented here is 
not perturbative. It consists in developing the stochastic variable as a 
power series in time, and using this to compute the short time expan- 
sion for the correlation functions. This, in turn, is extrapolated to large 
times and Fourier transformed to obtain the spectrum. A stochastic dia- 
gram technique is developed to facilitate computation of the coefficients 
of the relevant power series expansion. Two different types of long-time 
extrapolation technique, involving either simple exponentials or logarith- 
mic rational approximations, are evaluated for third-order diagrams. The 
analytical results thus obtained are compared with numerical simulations, 
together with exact results available in special cases. The agreement is 
found to be excellent up to and in the neighborhood of the critical point. 
The exponential extrapolation works especially well even above the criti- 
cal point at large N- values, where the dynamics is one of phase-diffusion 
in the presence of a spontaneously broken symmetry. A feature of this 
method is that it also enables the calculation of the steady state spectra 
of polynomial functions of the stochastic variable. In these cases, the final 
correlations can be non-bistable even above threshold, and the logarith- 
mic rational extrapolation has the greater accuracy. Finally, we emphasize 
that the technique is also applicable to more general stochastic problems 
involving spatial variation in addition to temporal variation. 

'Permanent address : School of Physics, University of Hyderabad, Hyderabad 500046 
(India) 
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1 Introduction 



Stochastic differential equations are a natural way of describing tfie interaction 
of a system witli a random reservoir. They were introduced by LangevinQ] to 
help explain Einstein's theory|^ of small particles immersed in a fluid, as ob- 
served by the biologist Robert Brown. More rigorous mathematical treatments 
were later introduced by Ito|^ and Stratonovich[Q. They have now diffused 
into many different areas Q of physics, chemistry and biology. In recent times, 
similar models have been utilized in ever more diverse fields, including engineer- 
ing, economics]^ , and even sociology. The essence of a stochastic differential 
equation is that it isolates a system of interest from the background of random 
events that may influence the system. Implicit in this formulation is the idea 
that the reservoir, or source of random fluctuations, evolves without reference to 
the system of interest. This simplifies the study of otherwise complex coupled 
phenomena. 

As an example, the calculation of correlation functions, and hence the spec- 
tra, of physical systems near phase transitions is of considerable theoretical inter- 
est. These have dynamical properties that are often conveniently described using 
stochastic differential equations. However, commonly used analytic techniques 
like linearization, frequently become invalid near phase-transition points. At 
the same time, while numerical simulation is possible, this is a time-consuming 
computational procedure without resulting in a great deal of insight. Thus, 
there is a need for techniques that give analytical results. Rather surprisingly, 
there are few systematic procedures for treating nonlinear stochastic differential 
equations under conditions where linearization is invalid. 

In this paper we consider the spectra of physical systems that are described 
by stochastic differential equations near a critical point phase transition. Such 
differential equations have a near universal applicability, both for equilibrium 
and non-equilibrium phase transitions For instance, they arise in theoret- 
ical treatments of single mode lasers ^ , inhomogeneously and homogeneously 
broadened two mode lasers ^ , and optical parametric amplifiers near thresh- 
old [0 A number of useful theoretical techniques jl^, ^ are known for 
these problems, some of which improve upon the the conventional Zwanzig-Mori 
projection operator method |p^ , p^ . However, these techniques are cumbersome 
for systems in higher phase space dimensions. 

Instead, we propose a simple, direct calculation which is based on the 
stochastic differential equation. The resulting integral expressions can be clas- 
sified diagrammatically, in a way that allows a straightforward calculation of 
essential combinatoric factors. The results give a power series in time which 
can be extrapolated to long times with reasonable accuracy in many cases. We 
analyze two possible extrapolation techniques, namely the exponential of a ra- 
tional function, and a series of simple exponential terms. Either method gives 
excellent results at or below the critical point. 

Above the critical point, we find differences in accuracy, and this can be 



2 



related to the dominant eigenvalue distributions for different types of equation 
and observable. Convergence is slowest when the spectrum has characteristic 
time-scales which are an exponential function of an equation parameter, as in 
the one-dimensional cubic stochastic process above threshold, which involves 
diffusion over a barrier. The exponential series method is preferable for simple 
types of spectra with only one or two dominant eigenvalues, which turns out to 
be the case for the iV-dimensional cubic stochastic process at large iV-values. 
The rational function method is best for complex spectra without much range 
in characteristic time-scales. An example of this is the intensity (x^) correlation 
spectrum of the one-dimensional cubic process, which can be represented with 
remarkable accuracy using low-order rational function extrapolation. 

More generally, we expect that this method can be applied to any stochastic 
differential equation where conventional linearization methods are inapplicable 
due to large nonlinear terms. Under these circumstances, it may be useful to 
have a nonlinear solution of the type derived here, as a starting point for a 
perturbative or asymptotic analysis. For these reasons, it is useful to analyze 
the simple cubic nonlinear case in detail, both as a test case for the stochastic 
diagram method, and as an elementary stochastic process of intrinsic interest. 

2 Stochastic Equations 

The method of stochastic diagrams to calculate solutions to stochastic differ- 
ential equations is normally applied in the frequency domain, where it corre- 
sponds to a perturbation theory expansion in a small coupling constant |20[ ^ . 
In these applications, there is a close resemblance to Feynman diagram tech- 
niques. In both cases, the starting point of the iterative method is the ap- 
proximate linearized solution to the problem, which becomes the solution to the 
entire correlation function in the limit as the coupling constant approaches zero. 
Frequency domain stochastic diagrams have many useful applications, includ- 
ing the stochastic quantization approach to quantum field theory. An essential 
difference between these methods and Feynman diagrams, is the appearance of 
stochastic terms that are averaged over at a later stage. 

We choose here to apply stochastic diagrams to the time domain correlation 
function. This has the advantage that there are no singularities in the series 
expansion coefficients, even at a critical point. A corresponding disadvantage is 
that the long time correlation functions cannot be directly calculated, and must 
be approximated by an extrapolation procedure that is based on some known 
property of the solution. In the examples given here, we use either simple expo- 
nentials or logarithmic rational function extrapolation, which results in analytic 
expressions for the approximate correlation function. A direct comparison with 
numerically calculated spectra will be used to demonstrate the great accuracy 
of this procedure in calculating spectra near critical points. It is less accurate 
above threshold in the bistable cases where stochastic 'barrier hopping' or 'tun- 
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neling' can occur, resulting in widely differing eigenvalues. This results in a 
reduced precision for the extrapolation. Methods based on multiple time scales 
can be used in these cases. 

Surprisingly, in the closely related higher-dimensional phase-diffusion prob- 
lem, stochastic diagrams give good results. An example of this is the laser above 
threshold. This is because the long time-scale here is not exponentially long, 
as it is in the tunneling cases. Good results are also found for non-bistable 
variables like the intensity, even when the underlying equations arc themselves 
bistable. This is because the extrapolation is carried out in terms of the cor- 
relation function, which has a different behavior to the underlying stochastic 
variable. In order to illustrate these various cases, we start with a very general 
form of stochastic differential equation. 

The equations we wish to treat are of the form 

x = A(x)+B(x).|(t) , (1) 

where the real noise sources S.i{t) have zero mean and are delta-correlated in 
time so that 

<e»(i)>=0 ;<C.m,{t')>^S,,Sit~t') . (2) 

Here x is a real vector of n-dimensions, A is an n-dimensional real polynomial 
vector function of x and B is an n x m dimensional real polynomial matrix 
function of x. The vector ^(t) is an m-dimensional real Gaussian stochastic 
process, interpreted in the Ito sense for computational simplicity. 

2.1 Iterative solutions 

The method of stochastic diagrams consists of performing an iterative solution 
for x(p) (t) so that 

x(rf(t)=xo+ rdt'[A(x("-i)(i'))+B(x("-i)(i'))-|(t')] , (3) 

Jto 

where x.'-°\t) = xq = x(to). Next, correlation functions of the typical form 

g^j{t,to) =< x,{t)xj{to) > - < x,{t) >< Xjito) > , (4) 

are evaluated to p*'' order, resulting in an expansion of Qij {t, tg) as a power series 
in T — t — to for t > 0. For any given term p in the power series, iterations 
must be carried out until all possible terms in are evaluated. The result 
of the iterations consists of integrals over time which will be represented as 
directed lines in the stochastic diagrams. In addition, there are polynomials in 
variables (represented as vertices), initial conditions in the variables (represented 
as terminating arrows, and treated as delta functions at the initial time) and 
stochastic terms (represented as crosses). 
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2.2 Cubic stochastic process 

Thus, for example, the solution of the well-known cubic stochastic process 

x^~x^+m , (5) 

has a first iteration, starting from a known initial value a;(0) — v at t = 0, of 

a;(i)(i) ^v + w{t)- [ dt'v^ . (6) 
Jo 

where w{t) = j^dt'^{t'). More generally, the n*'' iteration in this simple one- 
dimensional case is written 

x^''\t)=v + w{t)- f dt'{x'^''-^\t')f . (7) 
Jo 

Thus, we can expand the second term in the iteration as 

= y + yj(^t)- f dt'[v + wit')-t'vY 

Jo 

= v + w{t)+ f dt'[{t')\^-3it'fv^w{t')~3{t'fv^-6t'v^w{t') 
Jo 

+ 3t'v^w{t'f + 3t'v'^ -w{t'f ~3v^w{t')-3vw{t'f ~v^] . (8) 

We see that even this simple example leads to a large number of distinct terms, 
which need to be classified in a systematic way. In particular, while the leading 
term in the integral is of order t^, there are other terms of lower order present, 
including the stochastic terms and a term of order t which comes from the initial 
condition. 

2.3 Stochastic diagrams 

The next term in the iteration involves a cubic integral of x'--^\t), and clearly 
the combinatoric factors involved are more complex to three and higher orders. 
In order to simplify the counting of these factors, a diagrammatic classification 
can be introduced at this stage. In this classification, the terms are given di- 
agrammatically to first order in Fig (1). To second order, all possible terms 
in appear as 'legs' on the nonlinear vertex, to the next higher order, as 

shown in Fig (2). 

Not all terms will contribute to the same order in a power series in r, since 
the expectation value of a product of two stochastic integrals is proportional to 
r, while the product of two deterministic integrals is proportional to r^. This 
means that a reordering of the sequence is needed, to obtain a series of terms 
to a given order in r. The rules are simple: all vertices counts as one order in 
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r while stochastic terms count as half an order in r and initial values as zero 
order. If the terms in the reordered series are labeled as i^"^ (t) we can represent 
them according to Fig (3). On expanding all the relevant terms in Fig (3), we 
obtain: 

= V 

{t) = -^{t)+9v^ f wit')dt' + Qv^ f w{t')t'dt' 

Jo Jo 

= f «;2 it'Y'dt' + 9t;3 /* (t')dt' 
Jo Jo 

+ 18v^ f w{t')w{t')dt' - t^v'' - ^t^v'^ . (9) 
Jo 2 



^(5/2) 



'^{t)= f w'^{t')dt' . (10) 
Jo 



Here we have introduced the notation of: 

/■* 

lo 

Further rules in stochastic calculus (of the Ito variety) are that the expectation 
values of the products of initial terms with stochastic terms dccorrclatc to all 
orders at later times and all odd products of stochastic integrals average to zero. 
This means that the only surviving terms in the expectation value Qij must be 
the terms of integer order in the scries. For other types of expectation values 
(involving polynomials in x(t)), these extra terms must be retained. 

If we take expectation values of the relevant stochastic terms, they have the 
structure: 

{w\t)) = 

{wit)w(t)) = {w^t))=t^/2. (11) 

Hence, on decorrelating, integrating over time, and combining all the relevant 
terms up to third order we obtain (for the average and correlation function of 
x{t)): 

{x{t)) = {v-tv^ + ^t^{v^ -v) + h^^{llv^ -bv"^)) 

{x{t)x{Q)) = (?;2-to4 + ^t2(^6_^2)^ 1^3(^^^4_5^8)^^ (^2) 

These results are valid for an arbitrary initial distribution function. If carried 
out to higher orders, it is clear that they can describe either a transient process, 
or else a steady-state correlation in the time-domain, to any order in time. 
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3 N-dimensional cubic stochastic process 



Having introduced the stochastic diagram method, we now apply it to the N- 
dimensional cubic stochastic process 

Xi(t) = -rjiXi - fijkiXjXkXi + ^i{t) ; i,j,k,l = l,---,N . (13) 

Here summation over repeated indices is impHed. The coefficient fijki of the 
cubic terms can be taken to be symmetric in the the last three indices without 
any loss of generality. This stochastic equation, with appropriate choice of 
parameters accommodates the stochastic equations that have been considered 
in the context of single and two mode lasers and optical parametric amplifiers 
[1-6] . The quantities of interest are the equilibrium correlation functions 

4")(t) = lim [< x?(t + to)x^ito) >-< x^ir + to) X x]ito) >] • (14) 

This equation now has the added feature of a linear loss/gain term rj. When 
T] > 0, there is additional damping, and the system is below threshold in the 
usual sense. The deterministic critical point is at rj = Q. When rj < Q, the system 
has linear gain (like a laser above the lasing threshold), and the system is then 
above the critical point. However, it is worth noting that as the dimensionality 
increases, this type of classification which comes from a linearized analysis is 
rather misleading. In fact, the enlarged phase-space volume means that noise 
sources become increasingly important at large dimensionality - to the point that 
there is a reduced distinction between the above and below threshold cases. 



3.1 Steady-state behavior 

Steady-state behavior is most readily analyzed if, for simplicity, we confine 
ourselves to the following equation with A/'-dimensional rotational symmetry: 



x{t) = -r?x - x(x • x)/iV + ^{t) 



(15) 



This corresponds to defining the cubic coefficient as: fijki = [Sijdki + SikSji + 
Sii6jk]/iSN) . 

This stochastic equation has what is known as detailed balance - and hence 
an exact solution in the steady-state, found by examining the corresponding 
Fokker-Planck equation: 



d_ 

dt 



P{t,x)=LFpP{t,K)=Y^ 



d 

dxi 



^/N]xi 



1^ 
2d^ 



P{t,K) . (16) 



The equilibrium distribution is Pe(x) = Af ex.p[—V{x)], where V{x) is a 
potential function given by: 



y(x) = r?x • X (x • x.f/2N . 



(17) 
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The stochastic equation also has an exact relationship between the moments 
of different orders in the steady-state, which can be easily derived from the 
variable-change rules of Ito stochastic calculus. These are: 

^{2n + N- 2)M^"-^^ /{2N) - t]M^"') (18) 

Here we have defined A^*^"^ ([x-x]")e/iV, as a convenient normalized form 
of the moment. Although these recursion relations are useful, the mean-square 
fluctuations have to be calculated from the potential solutions given above. The 
quantity M^^'> = (x • x)e/7V can therefore be computed explicitly and is given 

m'^>^^u(^,^„)/u{^,Vn„) . (19) 

where U{a,x) denote the Whittaker functions For = 0, this expression 
simplifies to 

A<"'^,/|rm/rm . (20) 



3.2 = 1 case 

An important property of this potential in the one-dimensional case of = 1, 
is that it possesses a potential barrier at a; = 0, if 77 < 0. This means that there 
is a progression from a stable 'below-threshold' region for rj > 0, (where a; = 
is the deterministic stable point), to a critical region for 77 = characterized by 
large fiuctuations, and then to a bistable region for 77 < 0. This is characterized 
by local stability in the two potential wells at x = ±-\/|77|. 

The A^ = 1 case has been well-studied in terms of its eigenvalue spectrum. 
Any onc-dimcnsional Fokker-Planck equation can be transformed to an equiva- 
lent Schroedinger equation problem with imaginary times |ll^ , by introducing a 
Schrocdinger operator. In this case, it has the form: 

Ls = cxp[y(x)/2]LFpcxp[-y(x)/2] 
= [2V"{x)-{V\x)m+l^ 

= -i[x^^+,7x]2 + [r7 + 3a:2]/2+i^ . (21) 

At large positive values of 77, the corresponding Schroedinger potential re- 
duces to a harmonic oscillator problem, with quadratic potential. Transforming 
back to real time, the eigenvalues of the Fokker-Planck operator are of form: 



P„(x) = -A„P„(x) , (22) 



where A„ — nrj. Physically this is easy to understand. The equation is domi- 
nated by the linear decay rate 77, and integer multiples of 77 will occur through 
the decay of integer powers of the variable x. 
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At large negative values of rj, the equation is bistable, and there are two 
principle eigenvalues. A fast equilibration inside each potential well occurs, 
with an eigenvalue of A/ = 2|?7| in the limit of large t]. A slow decay also occurs 
through diffusion over the barrier. Ignoring terms in I77I of order (1) in the 
pre-factor, this gives a slow eigenvalue: 

ln[A,] ~ -Ay + ln[|r;|] ~ -r;V2 + ln[|r;|] . (23) 

It is significant for the stochastic diagram method, that this eigenvalue is 
exponentially small in the limit of large r]. Thus, we cannot expect an accurate 
estimate of the eigenvalue with any technique involving a finite scries of terms 
in 77, and any corresponding spectrum in which Ag is significant will not be able 
to be estimated with a finite expansion in powers of r]. 

3.3 Large N case 

For A'^ > 1, a similar progression from below to above the critical point holds 
deterministically, except that there is no bistable region. Instead, for 77 << 0, 
there is a region characterized by neutral stability in the subset of phase-space 
where |x| ~ y^\ri\N. Thus, there is a contimmm of possible deterministically 
stable behaviour. This phenomenon is sometimes called spontaneous symmetry 
breaking. To show this more clearly, consider the distribution Pij(i?), in the 
variable R = \x\'^ /N. This has a steady-state potential of Vz?.(i?) = N[riR + 
— (1/2 — 1/iV) ln(i?)]. As a result, for increasing N, the distribution in 
R is peaked more and more strongly near the value R^ = (i/2 + rf — ri)/2. In 
fact, due to the increase in phase-space volume as R increases, there is always 
an outward 'entropic' force even when -q > Q. This means that the stochastic 
equation at large N is not described well by the deterministic stability theory. 

In this limit, the radius approaches a fixed value, due to the balance between 
the outward entropic force due to increasing phase-space volume, and the inward 
force due to the nonlinearity. Thus, the moments A^^"' all factorise, and are 
given by: 

= ([v/2T^- 771/2)" . (24) 

The recursion relation for moments now simplifies, and it is straightforward to 
verify that the above solution does satisfy the recursion relation. 

Generally, in a stochastic equation, spontaneous symmetry breaking is ac- 
companied by a type of phase-diffusion, or tangential diffusion in a surface of 
dimension N — 1. Hence, the lack of bistability for any N greater than one 
resiilts in a dynamical behaviour in which diffusion still occurs, but with a re- 
duced dimensionality. These two types of above-threshold behaviour result in 
different dynamical regimes for the resulting correlation functions and spectra. 
In both cases, the above threshold dynamics typically involves more than dif- 
ferent time-scale. The radial relaxation to a stable point within a potential well 
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in the i?-space equations generally occurs much faster above threshold than the 
tangential diffusion. 

There arc corresponding changes in the large-N dynamics, and the physical 
explanation for this is interesting. In the limit of oo, the fast radial 

equilibration takes place in an approximately quadratic potential well at all 
values of 77. It may be noted that the corresponding (Ito) stochastic equation 
for radial equilibration involves noise in a multiplicative way. In the general 
case, we find that: 

R = -2r]R - 2i?2 + 1 + 2^^{t) , (25) 

where (C(i)^(i')) = S{t—t'). As well as having multiphcative noise, this equation 
also shows why the stochastic equation trajectories are confined to an increas- 
ingly small region in i?-space, as N increases. This occurs because the relative 
size of the noise term in the radial equation decreases as N increases. Thus, 
radial equilibration takes place with a fast relaxation rate of 

Xf = ARr, - 2r? = 2 V2 + rf . (26) 

In the tangential direction, diffusion takes place on a hyper-spherical surface 
of fixed radius defined by 

|x| = yiVi?; . (27) 

Suppose we define a coordinate system so the diffusion starts with a radial 
coordinate of xi = ^JNR^ at time t = 0. For small times, the other (tangential) 
coordinates obey the diffusion equation, so that 

{x-]) = t , (28) 

for J > 1. Since the radius is fixed by the radial equation, it follows that 
this corresponds to angular diffusion. Projecting each angular variable in turn 
onto a radius vector in a lower dimensional subspace reduces the length of 
the resulting vector. Finally, in the subspace of one dimension containing the 
original (starting) vector, we have: 

{xi{t)) ~ ^jNRn[l-t/NRr,](^-^) ~ ^V^exp [-(1 - l/N)t/2Rr,] . (29) 
This corresponds to a much slower tangential relaxation rate of 

A« = (1 - l/iV)/(2i?^) (30) 

in the large N limit. We note that this is not exponential in 77, unlike the one- 
dimensional case. Similar behaviour would occur in the case of finite N and 
large, negative r], which is also dominated by the tangential diffusion caused by 
spontaneous symmetry breaking. However, for the case of finite N values, the 
slow eigenvalue must reduce to Xg = r] in the limit of large enough positive rj. 
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4 N-dimensional stochastic diagrams 



The A''-dimensional equation clearly has the same structure as the integral equa- 
tion associated with the simple cubic process considered in the previous section, 
except for the linear terms which would complicate the diagrams if retained. In- 
stead, we can simply define yi{t) = exp{r]i{t — to))xi{t), which obeys a stochastic 
equation without a linear term. This can then be iterated as previously. The it- 
erative solution has the same diagrammatic structure as before. Using an initial 
estimate for Xi{t) of xf \t) — e~^^^'^~*°'>Vi, where the initial value is v = x(to), 
the basic iterative solution is given by; 

= e-'"(*-*°)x,(io) 

Jto 

(31) 

It is easily checked that replacing the approximations x^"', x^""^) by x, leads 
to an exact integral equation for x(t). 

We can now identify successive iterations with terms in the stochastic dia- 
grams for vector quantities Xi{t), where each vertex includes a term —fijkh and 
each directed arrow corresponds to / exp(— 77i(< — t'))... Thus, in evaluating the 
diagrams, each vector initial condition is replaced by Vj exp{—r]j{t — to)), and 
the noise term w{t) is replaced by: 

Wi{t)= f dtie-^'^*-''ki{ti) (32) 

. To order t^, this can be calculated using the diagrams in Fig (3a)- (3d), by 
making the associations given in Fig (4). Using the diagrams in Fig (3), one can 
easily derive an expansion for Xi{t) up to order t^. The details of the resulting 
stochastic integrals are straightforward, but rather lengthy. 

These results are given in the Appendix for the rotationally symmetric 
case. In the symmetric case, it is also clear that Q^^^^ (r) = if i ^ j and that 

= lim [< x(r-|-to)-x(fo) > - <x(T-hio) > • <x(io+T) >]/iV . 

(33) 

4.1 N-dimensional two-time correlation function 

Averaging the expression for Xi{t) thus obtained over ^(t)'s, expanding the ex- 
ponential factors and keeping all terms up to order , one obtains the result for 
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the steady-state or equilibrium two-time correlation function. In this expres- 
sion the two-time correlation is given in terms of the initial one-time moments 
of the stochastic process. Wc note that it is not essential, at this stage, to use 
equilibrium moments. The same general result occurs regardless of the initial 
condition, even for the case of transient correlations calculated without taking 
the steady-state limit. 

Wc will focus on the symmetric case here, and simplify the following expres- 
sion (derived in the Appendix), by using the definition of A^(") = < i?" >^ = 
< [x • x/N]"^ >e, where the subscript e denotes an equilibrium average: 



— T 



T 2N 



6 



6A^ 



+ 



13772 IN + 26 



6 



9?7 



(34) 



Next, we can substitute the known equilibrium moments to obtain a final 
expression for the correlation function in terms of the mean square fluctuation 
Al'^^ =< x-x/N >e, although still in a power series in r. Using the previous 
relations 7W("+i) = (2n + N- 2)7W("-i) /(27V) - r]M^''\ and defining 



a= 1/(27W«), 
one obtains the following power series: 



(35) 



— T 



l-ar + T'' 
N + 8 1 



127V-«+6^« 



N + 2 1 

AN ^ 2 

'4- 



77a 



127V 



It is convenient to re-express this as: 

3 

5(i)(r)=Al«[l-^a„T"] . 

where: 



(36) 



(37) 



n=l 



ai = a 

a2 = -{N + 2 + 2Nr]a) / (47V) 

as = ((7V + 8)a + 2Nri^a - (4 - N)r}) /(127V) 



(38) 
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Although the result is still expressed in terms of the correlation function through 
a, this quantity can be calculated exactly, either by integrating the distribution 
function numerically, or by using the Whittaker function representation. In the 
case of = 0, this reduces to a = ^yWJ8T {N/i)/T {{N + 2)/4). 

4.2 Correlations of polynomials 

The diagrammatic expression for Xi{t) in powers of r and W{t) can also be used 
to calculate the equilibrium correlations for any polynomial functions of Xi{t). 
Thus, in the case N = 1, for the equilibrium correlations for R = x^: 

g(2)e(T-) = < R(t + to)R{to) >- <R{t + to) >< R{to) > , (39) 

to— *oo 

we obtain 

^?(2)-(r) = [A^(2) _ (;^(i))2] _ ^[2A^(i)] + ^2j2 _ 2r,A^(i)] 

- r3[8A^W + ^r?'A^W] + 0(T*) . (40) 

Here we notice that A^^^) = x/2 — r]M.^^\ so the pre-factor in the above 
expression reduces to: 

[^(2) _ (^(i))2] ^ g(2)(o) ^ 2a^ -^ja - 1 _ (4^^ 

Just as in the case above, we can write the two-time correlation function 
down in terms of the individual power series terms, as: 

a(2) ^ 4a/[2a2_ 2770-1] 

a(2) ^ _ [8^2 _ 4^^] / [2a2 - 2r?a - 1] 

= 4a(4 + 2r?V3)/[2a^- 2770-1] (42) 

5 Spectral calculations 

The correlation function in the time domain must be extrapolated to long times 
in order to compute the spectrum, which involves a Fourier transform over all 
times. The general spectrum for any steady-state correlation function is: 

/•oo 

^("^(w) = 2i?e / dTg'^"\T)e"^^ . (43) 
Jo 

In order to perform the Fourier transform, some extrapolation of the power se- 
ries is required. In general, for an arbitrary initial condition, this is a difficult 
operation to perform. However, in the steady state, the un-subtracted correla- 
tion function must factorise at long times to the product of the mean values at 
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initial and final times. This means that the correlation function defined here 
gives rise to exponential decay at long times. 

In fact, for the type of stochastic differential equations considered here, we 
expect a discrete spectrum with an exponential decay at long times. However, 
a simple truncation of the power series at a finite order will not lead to an ex- 
ponential decay, so we cannot just truncate the power scries in time to obtain 
the spectrum. We will consider two different approaches to extrapolation. The 
first is to simply represent the correlation function with a finite series of ex- 
ponentially decaying terms, the second is to approximate the logarithm of the 
correlation function as a rational function. 

We assume that our starting point is an arbitrary correlation function Q, 
expressed as power series up to p-th order in the stochastic diagrams, of form: 

e(r)=c;(0)[l-^a„r"] . (44) 

71=1 

5.1 Simple exponential extrapolation 

This technique represents the correlation function as a finite series of decaying 
exponential terms. The coefficients can then be matched to the known power 
series in time on a term-by term basis. This method is especially useful when 
only a small number of eigenvalues dominates the spectrum. 

For a power series calculation to order r^, two distinct exponential terms are 
required. More generally, any correlation function expanded to order p = 2p' —1 
is represented using p' effective eigenvalues as: 

p' 

5(t) =e(0)^<?„exp(-A„T) . (45) 

71=1 

Here, for simplicity, we impose the restriction that J2n=i 9n — 1- 
obviously necessary that all the effective decay rates are positive. In the third 
order stochastic diagram case, on matching powers of r, one obtains: 



1 a\ + 3aia2 + Sas 

~ 2 2A 

1 a? + 3aia2 + 803 
92 = ;r + 



Ai = 



2 2A 
—A — 3a3 — a2ai 



1 



2a2 + a\ 



2a2 + al ^ ' 
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where the denominator term A is defined by: 



A = y Galas - 3af + 18010203 - 80^ + 9o§ . (47) 

5.1.1 Amplitude correlations - N=l case 

For the case of Q'^^^ at N = 1, the dependence of Ai and A2 on 77 is displayed 
in Fig (5a) by the two solid lines, where the upper line corresponds to Ai and 
the lower line to A2. There is a marked transition in this N = 1 case, between 
the non-bistable behaviour for >> 0, where the two time-scales are identical, 
and the bistable behaviour for ry << 0, where one time-scale becomes very long, 
corresponding to stochastic 'barrier-hopping' over the potential barrier in the 
distribution function at x = 0. In this region, the extrapolation technique used 
here is obviously less reliable, since the relevant eigenvalue is an exponential 
function of 77. One cannot accurately estimate these long-time tails on the 
correlation function, purely from the short-time information provided by the 
stochastic diagrams. In fact, there are other techniques based on multiple time- 
scales, which are more suitable in this above-threshold region. 

Nevertheless, the technique does generate the fast eigenvalue (A/ = 2|77|) 
correctly for large negative 77. For large positive 77, the harmonic oscillator 
predictions are regained. It is interesting to note that the fast eigenvalue in this 
case is A/ = S\ri\; this occurs because the symmetry of the problem means that 
even order eigenvalues are not significant in the dynamics of G^^-*(t) at large 
damping. The slow eigenvalue is not accurately reproduced at large negative 
?7, since this becomes exponentially slow (i.e, should be a straight-line graph). 
We will show later, by comparisons to numerical simulations, that the critical 
dynamics are reproduced accurately. 

5.1.2 Amplitude correlations - N=4 case 

In Fig (5b), the behaviour of the eigenvalues for Q^^^ at A'' = 4 is shown in 
the solid lines. Here we expect the slow eigenvalue to approach As = .375/(i?^), 
where i?,, = [■\/2 + Tf — T]\/2, while the fast eigenvalue should be A/ = 2-\/2 + rf. 
Since these are strictly large-N limits, we cannot expect these to be found ex- 
actly. These approximate results are actually reproduced with surprising pre- 
cision, especially in the phase-diffusion limit of large negative 77. Thus, we find 
As = .076 and A/ = 10.4 at = —5. The approximate predicted values would 
be As = .074, and \j = 10.6, with even better agreement at larger values of \rj\. 
At large positive ry, the slow eigenvalue approaches ry, and the fast eigenvalue 
approaches 3ry due to the term in the stochastic equation. This is a result 
due (as in the N = \ case) to symmetry properties; the intra-well eigenvalue 
with A = 2ry does not contribute to this correlation function. 
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5.1.3 Intensity correlations - N—1 case 



Finally, we consider a quadratic correlation function, which we term the inten- 
sity correlation. This is the case of G^^^ at N=l. These results for the relaxation 
rates arc given in the solid lines of Fig. (5c). Here wc sec no trace of the expo- 
nentially slow eigenvalue. This describes a sign-reversal process which has little 
or no effect on intensity correlations. Hence, all the observed relaxation rates 
arc caused by the higher-order eigenvalues for intra-wcU relaxation. Far below 
threshold, and above threshold, the eigenvalues approach 2\r]\ and 4:\r]\, which 
are characteristic of intra- well relaxation. Near the critical point of 77 = 0, there 
is strong critical slowing down, with longest time-scales (smallest eigenvalue) 
being found at 77 ~ —1.5. Although this region is bistable, it gives the slowest 
relaxation rate of As = 2.35; going further into the bistable region speeds up 
the intra-well relax;ation. 

5.2 Rational logarithmic extrapolation 

An alternative 'generic' technique, is to approximate the logarithm of the corre- 
lation function as a rational function, with a numerator of one order larger than 
the denominator. This is guaranteed to have an exponential behaviour at large 
T. We call this procedure rational logarithmic extrapolation. It can be applied 
to a power series of any order, as long as it is known that the series gives rise 
to exponential decay at long times. It is especially useful for complex spectra 
that may have several eigenvalues - as long as the values are not too different 
from each other. 

For a power series calculation to order r^, a quadratic rational function is 
required, so that we can approximate the correlation function to the given order 
as 



On matching powers of r with the previous power-series expression, one obtains 
the following general results: 



These can be used to obtain an extrapolated correlation function in any 
of the cases treated here. However, it is clearly important to ensure that the 
asymptotic coefficient, f3/j, is positive - otherwise no decay will occur. 




(48) 



a 



ai 

_2{a3 +aia2 + al/3) 

af -I- 202 
ai7 -I- al/2 + as 
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(49) 
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5.2.1 Amplitude correlations - N=l case 

In the expression for G^^\ the logarithmic expansion gives: 

ANa^ - QN-qa^ + 2Nrfa - 2{N - l)a - (4 - N)r] 
^ ~ 3N + 6 + GNrja - 6Na^ ' 

P ^ {SN + 6 + QN{r,a-a'))' 

7 12N{4Na^ - QNria? + 2N7fa - 2{N - l)a - (4 - N)r]) ' ^ ' 

The approximate expressions for the equilibrium correlations of x{t) given 
above arc characterized by two time scales - Xs = a — a and A; = 6/"f which 
govern the short and long time behaviors respectively. For the case of ^'^^ at 
N = 1, the dependence of As and A; on r] is displayed in Fig (5a) by the two 
dashed lines. As previously, there is a marked transition in this A = 1 case, 
between the non-bistable behaviour for 77 >> and the bistable behaviour for 
?7 << 0. Since the two time-scales here correspond to overall rates at long and 
short times, and not effective eigenvalues, the distribution of rates is different 
below threshold - eigenvalues with a low weighting do not contribute very much 
to the final rate. For this reason, we see no direct evidence for the fast time-scale 
below threshold, which corresponds to the relaxation of higher order eigcnstates. 

In the region of long time-scale 'barrier-hopping', the extrapolation tech- 
nique used here is less reliable. One cannot accurately estimate these long-time 
tails on the correlation function, purely from short-time information. This can 
be seen most clearly in the way that the slower of the two time-scales goes off 
the bottom of the logarithmic graph. At this stage, the longest time-scale is 
negative, indicating that the rational function approximation has broken down, 
and would predict an infinite or diverging spectrum. Obviously, the extrapo- 
lation is severely inaccurate at this point, and cannot be used this far above 
threshold. 

5.2.2 Amplitude correlations - N=4 case 

One might expect that the rational approximations used here should improve 
above threshold as N increases, as the equations are not bistable for large N. 
This hypothesis proves to be valid, as we show by the use of numerical stochas- 
tic techniques in the following sections. However, the improvement can already 
be seen in Fig (Sb). In Fig (5b), which gives Q^^^ at A = 4, the slowest relax- 
ation times above threshold are slightly too small compared to the exponential 
method, which indicates that the ratio of relaxation times is still too large for 
this method to give correct extrapolations, although the situation is much better 
than in the bistable case with A = 1. 
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5.2.3 Intensity correlations - N=l case 

For intrinsically non-bistable quantities like G^'^\t) the problems above do not 
occur at all, as shown in Fig (5c), which graphs G^'^^ at N = 1. Here the results 
of both extrapolation methods give similar behaviour. This is not immediately 
evident from the graphs, as the rates defined here do not have identical interpre- 
tations. In fact, the generic method of rational function extrapolation is actually 
better than the exponential method in this case. In order to demonstrate this 
in detail, we must turn to the full spectral calculation, which include both the 
relaxation rates and the relative weights. 

Before turning to the spectral results, we note that the two-time correlation 
function for the quadratic correlations in the rational approximation have quite a 
simple form. Extrapolating to large r using the rational function approximation 
as above, we obtain 



((1)121 



exp - 



1 + c'r 



where, using the result Ad^"^^ = 1/2 — r]A4^^\ 
a' = 4A^(i)/(l-2j7X(^)-2[>t(i)]2) 
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(51) 



rjM 



(1) 



(52) 



It is interesting to note from Fig(5c), that even for N = 1, the variable 
is clearly not bistable, and shows no sign of the characteristic long time-scales 
of bistable variables. Instead, there is a critical slowing-down near r] = 0, with 
shorter time-scales at all other r] values. This also agrees with the exponential 
extrapolation results. 



5.3 Spectral results 

Having computed a long-time extrapolation to the two-time correlation function, 
it is now possible to calculate the spectrum. A simple Lorentzian spectrum is 
generated by the first order stochastic diagrams. This approximation we shall 
see is surprisingly close to the correct spectriim even for finite N values, es- 
pecially at high frequencies. The reason for this is that the middle to high 
frequency spectral behaviour is mostly due to the change in slope in the two 
time correlation function near r = 0, due to the fact that the steady state corre- 
lation function must be a function of |r|. The low frequency spectral behaviour 
near a; = has additional contributions due to the r — > oo behaviour of the 
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correlation function, which is not always given accurately using a first order 
expansion in r. The lowest order spectral contribution is therefore 



S'-^\u) = 2g(i)(0)i?e / dTe~'^l"l+*'"" = ^ , (53) 

Jo a + ^ 

where a = 1/[2{R)] is defined as in the previous sections. Thus, in the large N 
limit we expect to find that: 

lim = 7 , " . (54) 

N^oo ^ ' (l + (2i?,,w)2 ^ ^ 

To higher order, it is necessary to choose which extrapolation method to use. 
The procedure of extrapolating the correlation function to large t by expanding 
it as a series of exponentially decaying terms, is simple to Fourier transform, 
since this clearly results in: 

'^M-^wf:!^ . (55) 

This results in a Fourier transform as a sum of Lorentzian components. For the 
three previous cases treated of = 1, = 4, and the intensity spectrum for 
= 1, the spectrum thus obtained is plotted in Fig (6a) -(6c) for various values 
of 77. The clear progression from bistable behaviour, to spontaneous symmetry- 
breaking, to non-bistable behaviour, is shown in these three graphs; as the peak 
spectral intensity is greatly reduced for the spectra with shorter characteristic 
time-scales. 

While this procedure has advantages as far as calculating the Fourier trans- 
form is concerned, it is not always the best extrapolation. Accordingly, we also 
consider a straightforward numerical Fourier transform, which can be used to 
calculate the spectrum from the rational function approximation to the corre- 
lation functions. 



5.4 Direct numerical simulation 

As there are no exactly known analytic expressions for these spectra, we have to 
resort to stochastic numerical techniques to check the accuracy of the spectrum 
using the truncated diagram method. Thus, in order to determine the accuracy 
of the correlation function, a direct simulation of the relevant stochastic differ- 
ential equations was used. These simulations employed a strongly convergent 
semi-implicit numerical algorithm , with checks on both numerical sampling 
error and truncation error due to finite step-size. The actual algorithm used 
employed four iterations of the nonlinear implicit equations at each step. After 
starting the trajectories in a Gaussian distribution with variance equal to the 
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known steady state variance, each trajectory was integrated for a total elapsed 
time much longer than the correlation time. 

Typically, t = 100 — 400 was the maximum time used, with longer times 
being employed for calculations above threshold, where the correlation time 
is longer. In the simulations, the total number of trajectories employed was 
10®, in order to reduce the sampling error to typically about 0.1 — 0.2% at 
low frequencies - although this typically rose to about 0.5% above threshold, 
presumably because of the highly non-Gaussian individual trajectory statistics 
in these cases. Sampling errors were checked by subdividing the results into 
1000 sub-ensembles, which were numerically Fourier transformed and averaged 
individually. The spectral results were then averaged over the sub-ensembles, 
and the error-bar in the overall mean was estimated using standard Gaussian 
distribution error estimates - on the basis that sub-ensemble means should have 
a Gaussian distribution according to the central limit theorem. 

By reducing the time step to a small value (typically At = 0.01 — 0.05), 
the errors due to the finite step-size were typically kept to below 0.5%. This 
was estimated by calculating all spectra at two different step-sizes, but with the 
same underlying stochastic noise terms, and comparing the results. The two 
error-bars were added to give the final numerical error-estimates. No significant 
error from finite spectral windowing was found, although this would be expected 
to give problems in the extreme bistable regions. 

By comparing with exact zero-frequency results, this technique of error- 
estimation proved a reliable method, in the sense that the discrepancies were of 
the expected size. 

Thus, for example, the numerical spectrum for the case of = 1 at threshold 
(77 = 0) has the calculated value of S{0) = 0.966 ± 0.007, using 10® trajectories, 
a window of tMAX = 200, and a step-size of At = 0.025. The error due to the 
finite step-size contributes ±0.005 to the total error. The corresponding exact 
result, as explained in the next section, is S{0) = 0.975 - giving a slightly greater 
actiial numerical error than tlic^ estimated one standard-deviation error-bar. By 
comparison, the extrapolated exponential series analytic result is S{0) = 0.970, 
which is very close to the exact resiilt. 

For the case of A'^ = 1 above threshold (77 = —1.5) the calculated numerical 
simulation value is S{0) = 10.01 ± 0.10, using 10® trajectories, a window of 
tMAX — 400, and a step-size of At = 0.05. Because the step-size is relatively 
large (in order to maximize the time- window) , the error-bar in this case is mostly 
due to the finite step-size, which contributes an error of ±0.08 to the total error. 
The corresponding exact result is S{0) = 10.11 - within the estimated error-bars. 
By comparison, the extrapolated analytic result is S(0) = 9.06, which gives an 
increased extrapolation error, as expected. 

The resulting numerical estimate of 5^"^ (w) was: 

Si:Li^) = lim < > /T , (56) 

1 — >oo 
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where the Fourier transform Aa;"(a;) is defined as 

A5r(c^) = r die-*[^rw - i^nm ■ (57) 

5.5 Compcirison of results 

The results of the two procedures, i.e. stochastic diagrams and numerical sim- 
ulations, are compared in Figs (7), (8) and (9) for iV = 1, iV = 4 and the 
intensity spectrum with N = 1 respectively. In each figure there are four lines, 
which arc the Lorentzian approximation (dotted line), the rational approxi- 
mation (dashed-dotted line), the second order exponential series extrapolation 
(dashed line) and the direct numerical simulations (solid line). 

In each figure there arc sub-figures which correspond to different values of 
the driving field rj, which are taken through a range of values from far below 
threshold to above threshold. We notice that agreement is generally excellent 
(close to the simulation error-bars) for all methods below threshold. Errors are 
always worst at low frequencies, where the results are the most sensitive to the 
extrapolation error at long times. They are also worst for the single exponential 
extrapolation than the higher-order extrapolation methods, as expected, and 
usually best for the exponential series method. 

5.5.1 Amplitude correlations - N=l case 

Below threshold, Fig (7a) shows the four different spectral results near zero- 
frequency, thus giving a magnified view of the results. Errors are much smaller 
at higher frequencies, where all the techniques agree to within the simulation 
error-bars. It can be seen that the exponential series method (dashed line) gives 
the best low- frequency agreement to the simulation (solid line). The residual 
difference is about equal to the intrinsic numerical discretization and sampling 
errors, while the other two methods give small, but marginally significant dis- 
crepancies. 

In Fig (7b), at the critical point of = 0, the exponential series gives a 
prediction at iV = 1 of S{0) = 0.9702. This is in quite remarkable agreement 
with the simulated value of S{0) = 0.966 ± 0.007. By comparison, the other two 
methods are again either significantly higher (rational logarithmic), or lower 
(Lorentzian approximation). 

Above threshold, however, in Fig (7c), we see that the agreement is outside 
the error-bars even for the exponential series method, thus indicating a reduced 
accuracy in the long-time extrapolation. This is expected, in view of the fact 
that the long-time eigenvalue is an exponential function of 77 - rather than a finite 
algebraic expression, as would be generated from the stochastic diagrams. Here 
the errors increase to about 15% for r] = —1.5, at zero frequency, in the rational 
approximation, with much worse errors in the single exponential approximation. 
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However, in this case the exponential series still gives the best result, with an 
error of less than 10% . 

5.5.2 Amplitude correlations - N=4 case 

Fig (8) shows the spectrum of x{t) as a function of 77 at = 4. In Fig (8a), 
with r] = 1 (below threshold) the difference between the various approximations 
and the numerical simulation is about equal to the sampling error-bars obtained 
with 10^ trajectories. Thus, it is not possible to distinguish between any of the 
extrapolations in this case. At threshold, in Fig (8b), the differences increase, 
and the exponential series method is clearly better than a single exponential 
extrapolation. Above threshold, in Fig (8c), the accuracy of both the single ex- 
ponential and the rational extrapolation diminishes, relative to the exponential 
series method - due to the two dominant eigenvalues in this case. For A'^ = 4, 
the above threshold error with the rational extrapolation reduced to 7%, as the 
multiple time-scale problem is less significant in this case. The error in the ex- 
ponential series method is less than 1%, i.e., of the same order as the intrinsic 
sampling errors. 

5.5.3 Intensity spectrum - N=l case 

Fig (9) shows the spectrum of x'^{t) as a function of rj as given by the dif- 
ferent approximations. We see that there are obvious differences between this 
and the previous cases. As x^{t) is not bistable, the spectrum does not have a 
large 'spike' as 77 ^ —00. This means that, unlike the previous examples, the 
agreement between the rational function extrapolation and simulation methods 
remains of the order of the numerical sampling error-bars even above thresh- 
old. Since it would presumably require more than 10^ trajectories to resolve 
the differences in these spectral results, we have not attempted to accurately 
estimate the extrapolation errors here. For this problem, it is clear that the 
analytical theory is rather competitive with numerical simulations, which are 
always subject to numerical sampling error. The nearly perfect agreement for 
rational function extrapolation is a surprising result, given that it is obtained 
from only a small number of stochastic diagrams. In this case, the exponential 
series method gives slightly worse results, presumably because there are several 
competing eigenvalues rather than just two. 

5.6 Exact results 

Although the agreement between the analytic results and the numerical simula- 
tions is generally excellent (except for bistable variables) there are some features 
worth a closer examination. Firstly, we emphasize again that discrepancies are 
only evident at low frequencies. This is simply because the correlation function 
tails, although only contributing a small part of the spectrum, cannot always 
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be accurately extrapolated from the short time power series expansion. This 
discrepancy at or below the critical point is always small (of order of 1 — 2% at 
the critical point rj = 0). The error is naturally larger in the simpler Lorentzian 
approximation, but in no case does it exceed 5% with i] = 0. Secondly, the low 
frequency extrapolation error is larger above threshold as rj — > — oo, especially 
in the bistable case of x{t) spectrum with = 1. This result is expected, since 
calculation of the slow eigenvalue in this case requires an infinite series in rj. 

In the particular case of = 1, the zero frequency spectrum is known 
exactly [7] and is given by 

S(-\0)^A r dxl^^p^ , (58) 

where 

/("'(.t) = - r dxiK- < >]Pe,(xi) . (59) 

J — oo 

Here Peq{x) denotes the normalized equilibrium distribution. The discrepancy 
between this exact result at = 1, and the approximate spectra, is given in 
Fig (10). 

In the bistable case of S'^^\ it is clear that the error increases very rapidly 
as 77 — *■ — cxD, especially for the single exponential and rational approximations. 
The reason for this is due to the well known bifurcation of the deterministic 
steady state in this case, when 77 << 0. The system can only switch from one 
solution near a; = 77 to the other near x = on exponentially long time scales. 
This can be seen on comparing Fig (11), which is computed at the deterministic 
threshold of 77 = 0, with an above threshold numerical simulation in Fig (12). 
In this region and above, a multiple time scale technique would give the best 
results, with other methods being used to estimate the slow eigenvalue. As we 
are interested here in the critical region around 77 = 0, we simply note this 
problem here, rather than pursuing it in further detail. Techniques for dealing 
with multiple time-scales of these barrier-hopping equations were first developed 
by Kramers p5| , and are quite well understood. 

This type of problem is greatly reduced for A^ > 1, where there is a tan- 
gential diffusion above threshold, rather than barrier-hopping. This results in 
an increasing accuracy of the present technique for large A^, as the system dy- 
namics reduces to just two time-scales, both of which have finite expression 
in terms of 77. In this limit, the behaviour is closely analogous to the well- 
known problem of diffusion in a curved space. While there are no examples of 
physical systems with such large-dimensional order-parameters, many physical 
systems (like the laser above threshold, or Bose-Einstein condensates) display 
this type of spontaneous symmetry breaking accompanied by tangential diffu- 
sion in phase-space. The case chosen (A^ = 4) is typified by a two mode laser 
with inhomogeneous broadening, so there is no strong mode-competition. As 
we have seen, the stochastic diagram method gives very accurate results when 
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compared with numerical simulations. However, we have no exact result in this 
case. 

Multiple time-scale problems do not occur in the case of non-bistable quan- 
tities like x^{t) whose spectrum is given in Fig. (10b). The exponential series 
method gives good agreement, as expected. Agreement with the rational func- 
tion extrapolation in this case is excellent even above threshold. It is so good 
that the exact spectral result at w = cannot be told apart from the approxi- 
mate value, so we have not included the rational approximation in this graph. 

6 Summary 

The stochastic diagram technique is a method of classifying iterative terms in 

a stochastic power scries expansion in time, so that terms of the same order 
in time are grouped together. This involves analyzing the deterministic and 
stochastic order, since they give rise to different types of contributions. 

In this paper we have used this technique to analyze critical amplitude and 
intensity spectra, by considering either exponential series or logarithmic ratio- 
nal function extrapolation. The results are inherently non-perturbative, and are 
very accurate except in bistable regions. Thus, the results arc valid at the criti- 
cal point, where linearized spectra diverge. Good results are also obtained even 
above threshold using exponential series extrapolation, when there is sponta- 
neous symmetry-breaking. The technique works especially well for intrinsically 
non-bistable quantities like the intensity or radial spectrum. In these cases, we 
find that rational function extrapolation is of greater accuracy. 

The general feature of these results for the case of the cubic stochastic pro- 
cess, is that we see the expected critical slowing down at the threshold or crit- 
ical point at ?7 = 0. Below threshold, there is rather stable behaviour, except 
at high dimensionality where the outward entropic effects of the phase-space 
volume are increasingly important. Above threshold, there is a bistable region 
for A'' = 1, where our methods do not converge quickly, due to the exponentially 
long time-scales involved. For this region, asymptotic methods involving mul- 
tiple time-scales would be more suitable. For A^ > 1, our method gives much 
better results than for A' = 1 above threshold, since the dynamics in this region 
are then dominated by a type of phase- diffusion at fixed radius, rather than by 
barrier-hopping. 

The problems treated here were of a relatively simple type, to allow com- 
parisons with exactly known results. We emphasize that the stochastic diagram 
approach is rather general. It is not restricted to equilibrium correlations, al- 
though extrapolation is not always possible in transient calculations. Nor is 
it restricted to problems just involving temporal variation: partial stochastic 
differential equations can also be treated in this way. 

In summary, this novel stochastic diagram technique appears to have many 
possible applications for calculations involving stochastic differential equations. 
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Appendix 



In this appendix we briefly outline the details of the calculations leading to the 
result given in (34), using the rules given in Fig (4). 

The contributions to the diagrams in Fig (3a) can easily be written down as 
follows 



Si^^'^\t) = Wi{t) 
Next, evaluating the diagrams with one vertex in Fig (3b), gives: 
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The next order terms in Fig (3c) include both one- vertex terms with additional 
noise factors, and two-vertex nested integral terms: 
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Finally, the relevant third order terms are: 
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In the above equations summation over repeated indices is implied. The next 

step consists in 

(a) multiplying the above expressions by Vi = Xi(0) and summing over i 

(b) averaging the resulting expressions over the noise sources and the initial 
values and adding up all the contributions. 

This leads to 
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(5) 



For the symmetric case with fi^u given by fijki = [5ij5ki+5ik5ji+5ii5jk]/ {iN) , 
the above results can be simplified. This is not essential to the method, and 
neither is the use of a steady-state initial condition at this stage. On explicitly 
carrying out the summations, we obtain 
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(6) 



Finally, replacing the initial averages by equilibrium averages and expanding 

the time integrals in powers of t we obtain 



(xW-x(O)) 
N 



+ 
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TV 

(iV + 2) 
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TV 



(7) 



On collecting coefficients of like powers of t we are led to the expression (34) 
for the short time expansion of the two-time correlation function in the N- 
dimensional case with rotational symmetry. This also agrees with the earlier 
result of Eq (12), for the simpler case of A'' = 1 and ry = - as expected. 
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Figure 1: Diagrammatic representation of the iterative solution of the cubic 
stochastic equation to zeroth and first order. 

Figure 2: Diagrammatic representation of the iterative solution of the cubic 
stochastic equation to second order. 

Figure 3: Diagrammatic representation of the iterative sohition of the cubic 
stochastic equation reordered by collecting together the diagrams which con- 
tribute to the same power of r. The power of r that each diagram contributes, 
is equal to the number of vertices plus half the number of crosses. Figs (a) 
-(d) indicate successively higher order diagrams, including fractional stochastic 
orders. 

Figure 4: The rules for writing down the contribution of a given diagram. 

Figure 5: The dependence of the apparent relaxation rates on the loss rate 

rj as given by the exponential series (solid line) and rational function (dotted 
line) extrapolations. Note that the two types of extrapolation implicitly define 
different types of relaxation rate. In each graph, the lower of the two lines 
corresponds to the slower relaxation rate, the upper to the faster relaxation 
rate. Cases treated are for (a) G^^\ N = 1; (b) G'^^\ N = 4; (c) G^^), N = 1. 

Figure 6: The approximate equilibrium spectrum for range of values of r], using 
the exponential series method. Cases treated are for (a) S^^\u)), N = 1; (b) 
5(i)(w), iV = 4; (c) 5(2)(w), iV= 1. 

Figure 7: Comparison of the equilibrium spectrum S'^^^uj) of x{t) for N = 1 
as given by stochastic numerical simulations (solid line), exponential series ex- 
trapolation (dashed line), rational function extrapolation (dashed-dotted line), 
Lorentzian approximation (dotted line) for (a) = 1, (b) 77 = 0, (c) rj = —1.5. 
Error-bars for discretization and sampling errors at zero frequency are: (a) 
AS = 0.002 , (b) AS = 0.007 , (c) AS = 0.1. 

Figure 8: Comparison of the equilibrium spectrum S^-^^oj) of x{t) for iV = 4 
as given by stochastic numerical simulations (solid line), exponential series ex- 
trapolation (dashed line), rational fimction extrapolation (dashed-dotted line), 
Lorentzian approximation (dotted line) for (a) = 1, (b) r; = 0, (c) rj = —1.5. 
Sampling error-bars at zero frequency are: (a) AS = 0.002 , (b) AS = 0.01 , 
(c) AS = 0.1. 
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Figure 9: Comparison of the equilibrium spectrum S^'^^ui) of x'^{t) for 
A'^ = 1 as given by exact numerical simulations (solid line), exponential se- 
ries extrapolation (dashed line), rational function extrapolation (dashed-dotted 
line), Lorentzian approximation (dotted line) for (a) 77 — 1, (b) 77 = 0, (c) 
rj = —1.5. Sampling error-bars at zero frequency are: (a) AS^^) = 0.0003 , (b) 
A5(2) = 0.001 , (c) AS'^^) = 0.004. 



Figure 10: Comparison of the equilibrium spectrum for = 1 at zero fre- 
quency as given by the exact analytic expression (solid line), exponential se- 
ries extrapolation (dashed line), rational function extrapolation (dashed-dotted 
line), Lorentzian approximation (dotted line) for a range of values of 77. Cases 
treated are: (a) 5(i)(0); (b) S<^^\0). 



Figure 11: A sample numerical run showing critical fluctuations of x{t) at 
threshold, for the case r] = 0, N = 1 . 



Figure 12: A sample numerical run showing bistability of x{t) above threshold, 
for the case r] = —1.5, N = 1 . 
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Fig. 1 
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X(l/2)(t) 

Fig. 3a 



Fig. 3b 
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Fig. 3c 
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Fig. 3d 
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t = -f^kijW exp[-ri.(t-t')] [ . . . ] 



Fig. 4 
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